Lab 8: Multilevel Models

Variation explained

Let’s first run a null multilevel model. Doing so will give us estimates of the amount and percent of variation explained at each level - tract and metropolitan area. To run a multilevel model, we’ll use the lmer() function, which is a part of the lme4 package.

library(lme4)

Next, run an empty random intercept model, which does not include any tract or metropolitan area level characteristics.

re.int.model1 <- lmer(lhinc~ (1|GEOIDm), data=ca.metro.tracts)
summary(re.int.model1, correlation = FALSE)

## Linear mixed model fit by REML ['lmerMod']
## Formula: lhinc ~ (1 | GEOIDm)
##    Data: ca.metro.tracts
## 
## REML criterion at convergence: 9368
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.6789 -0.6673 -0.0270  0.5731  3.3935 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  GEOIDm   (Intercept) 0.03683  0.1919  
##  Residual             0.19124  0.4373  
## Number of obs: 7841, groups:  GEOIDm, 26
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 10.85216    0.03892   278.8

I use correlation = FALSE in the summary() function to suppress the Correlation of Fixed Effects results, which frankly I have a hard time interpreting. The lmer function’s arguments take on a similar form as lm() except the use of (1|GEOIDm). The argument (1|GEOIDm) indicates the random effects of the equation, including the second level (metropolitan area) variation. The following table shows how to specify the different multilevel models using this ( | ) format.

There are other specifications that deal with different types of nesting structures. You’ll explore one of these alternative nesting structures in homework 4.

Going back to our output of results, the results under the “Random effects” heading provide the variation at each level - level 1 tract (Residual) and level 2 metro area (GEOIDm). We actually want the percent of variation explained, so we have to divide the RE of GEOIDm by the total variation (Residual + GEOIDm). We can do this by using the re_var() function, which extracts the random effects of a saved multilevel object. This function is in the sjstats package, which we’ll need to load.

library(sjstats)
var.est<-re_var(re.int.model1)
var.est[2]/sum(var.est)

## GEOIDm_tau.00 
##     0.1614856

The proportion explained at each level is known as the intra-class correlation (ICC). The icc() function in the sjstats package calculates the ICC directly for the higher level(s).

icc(re.int.model1)

## 
## Intraclass Correlation Coefficient for Linear mixed model
## 
## Family : gaussian (identity)
## Formula: lhinc ~ (1 | GEOIDm)
## 
##   ICC (GEOIDm): 0.1615

We find that the proportion of variation in median Hispanic Household income explained by differences across metro areas is 16%. The rest is explained by census tracts and potentially other levels.

Note that our outcome lhinc is continuous. If we have a non-continuous outcome such as a count or categorical variable, we use the glmer() function, and specify the outcome distribution using the family argument. For example, if you have a binary outcome, specify family = binomial, which runs a multilevel logit regression model.

Random intercepts with covariates

An empty random intercept model gives us percent variation explained at each level. Let’s now add covariates to try to find which specific variables might explain variation in the outcome. First, let’s add the tract-level variables

re.int.model2 <- lmer(lhinc~ concd  +pfb + phhown   +phmov + phisp 
              +  (1|GEOIDm), data=ca.metro.tracts)
summary(re.int.model2, correlation = FALSE)

## Linear mixed model fit by REML ['lmerMod']
## Formula: lhinc ~ concd + pfb + phhown + phmov + phisp + (1 | GEOIDm)
##    Data: ca.metro.tracts
## 
## REML criterion at convergence: 3839.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -10.9783  -0.5440   0.0242   0.5689   4.4519 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  GEOIDm   (Intercept) 0.01994  0.1412  
##  Residual             0.09414  0.3068  
## Number of obs: 7841, groups:  GEOIDm, 26
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept) 10.695049   0.032077 333.420
## concd       -0.200106   0.006615 -30.250
## pfb         -0.198233   0.033879  -5.851
## phhown       0.712921   0.015280  46.657
## phmov        0.032258   0.033087   0.975
## phisp       -0.312943   0.020315 -15.405

If you are annoyed by the results shown in scientific notation, use the scipen=alpha argument in the options() command, where alpha is the maximum number of digits for the result to be still expressed in fixed notation.

options(scipen=1)
summary(re.int.model2, correlation = FALSE)

## Linear mixed model fit by REML ['lmerMod']
## Formula: lhinc ~ concd + pfb + phhown + phmov + phisp + (1 | GEOIDm)
##    Data: ca.metro.tracts
## 
## REML criterion at convergence: 3839.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -10.9783  -0.5440   0.0242   0.5689   4.4519 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  GEOIDm   (Intercept) 0.01994  0.1412  
##  Residual             0.09414  0.3068  
## Number of obs: 7841, groups:  GEOIDm, 26
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept) 10.695049   0.032077 333.420
## concd       -0.200106   0.006615 -30.250
## pfb         -0.198233   0.033879  -5.851
## phhown       0.712921   0.015280  46.657
## phmov        0.032258   0.033087   0.975
## phisp       -0.312943   0.020315 -15.405

Now include the metropolitan area characteristics.

re.int.model3 <- lmer(lhinc~ concd  +pfb + phhown   +phmov + phisp 
              +lhincm + punempm + hwdis + fbdis + (1|GEOIDm), data=ca.metro.tracts)
summary(re.int.model3, correlation = FALSE)

## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## lhinc ~ concd + pfb + phhown + phmov + phisp + lhincm + punempm +  
##     hwdis + fbdis + (1 | GEOIDm)
##    Data: ca.metro.tracts
## 
## REML criterion at convergence: 3793.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -10.9789  -0.5422   0.0271   0.5666   4.4459 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  GEOIDm   (Intercept) 0.00187  0.04324 
##  Residual             0.09414  0.30683 
## Number of obs: 7841, groups:  GEOIDm, 26
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  4.120523   1.058315   3.893
## concd       -0.198102   0.006595 -30.040
## pfb         -0.199801   0.033701  -5.929
## phhown       0.712823   0.015254  46.729
## phmov        0.034829   0.033045   1.054
## phisp       -0.317133   0.020129 -15.755
## lhincm       0.597442   0.090403   6.609
## punempm      0.951119   0.623915   1.524
## hwdis        0.821212   0.131610   6.240
## fbdis       -1.104715   0.479702  -2.303

The tract-level coefficients don’t change much in terms of their effect sizes and statistical significance compared to the single and fixed effects regression results. However, we do see significant changes in the metro-level characteristics compared to the single-level results. The direction of the effects do not change, but foreign-born segregation is now significant at the 5 percent level. Moreover, the effect sizes are much larger in the multilevel model. How would you explain these changes?

Random coefficients

We might be interested in determining whether the effects of certain tract-level variables vary across metropolitan area. Here, we run a random coefficient or slope model. First, run a random coefficient model for phhown without a random intercept.

re.slope.model1 <- lmer(lhinc~ concd  +pfb + phhown   +phmov  + phisp 
              +lhincm + punempm + hwdis + fbdis + (0+phhown|GEOIDm), data=ca.metro.tracts)
summary(re.slope.model1, correlation = FALSE)

## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## lhinc ~ concd + pfb + phhown + phmov + phisp + lhincm + punempm +  
##     hwdis + fbdis + (0 + phhown | GEOIDm)
##    Data: ca.metro.tracts
## 
## REML criterion at convergence: 3810.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -10.9936  -0.5446   0.0240   0.5680   4.4731 
## 
## Random effects:
##  Groups   Name   Variance Std.Dev.
##  GEOIDm   phhown 0.01032  0.1016  
##  Residual        0.09424  0.3070  
## Number of obs: 7841, groups:  GEOIDm, 26
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  2.024789   0.797788   2.538
## concd       -0.196485   0.006605 -29.746
## pfb         -0.174621   0.033512  -5.211
## phhown       0.712717   0.030367  23.470
## phmov        0.030580   0.033068   0.925
## phisp       -0.326575   0.019913 -16.400
## lhincm       0.775635   0.067472  11.496
## punempm      0.862818   0.546285   1.579
## hwdis        0.597971   0.072915   8.201
## fbdis       -0.003154   0.341459  -0.009

The results under the Random effects section tells us that the random slope for phhown has a variance of 0.01032 - that is, the amount of variability in the slope of phhown across metro areas is 0.01032. How do we know whether this variability is statistically different from 0? We can use the function ranova() in the lmerTest package, which uses a Chi-square test to test against the null of variance equal to 0. First, load the package

library(lmerTest)

Plug in your saved random slope model above into ranova()

ranova(re.slope.model1)

## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## lhinc ~ concd + pfb + phhown + phmov + phisp + lhincm + punempm + 
##     hwdis + fbdis + (0 + phhown | GEOIDm)
##                                 npar  logLik    AIC     LRT Df Pr(>Chisq)
## <none>                            12 -1905.4 3834.8                      
## phhown in (0 + phhown | GEOIDm)   12 -1896.5 3817.1 -17.727  0

We can add a random intercept by replacing the 0 in ( | ) by 1.

re.int.slope.model <- lmer(lhinc~ concd  +pfb + phhown   +phmov  + phisp 
              +lhincm + punempm + hwdis + fbdis + (1+phhown+pfb|GEOIDm), data=ca.metro.tracts)
summary(re.int.slope.model, correlation = FALSE)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## lhinc ~ concd + pfb + phhown + phmov + phisp + lhincm + punempm +  
##     hwdis + fbdis + (1 + phhown + pfb | GEOIDm)
##    Data: ca.metro.tracts
## 
## REML criterion at convergence: 3687.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -11.0409  -0.5312   0.0189   0.5638   4.5105 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr       
##  GEOIDm   (Intercept) 0.02643  0.1626              
##           phhown      0.01365  0.1168   -0.74      
##           pfb         0.13007  0.3606   -0.95  0.49
##  Residual             0.09244  0.3040              
## Number of obs: 7841, groups:  GEOIDm, 26
## 
## Fixed effects:
##                Estimate  Std. Error          df t value Pr(>|t|)    
## (Intercept)    4.865846    0.763619   81.047043   6.372 1.06e-08 ***
## concd         -0.199236    0.006601 7278.598287 -30.182  < 2e-16 ***
## pfb            0.164583    0.091567   23.694727   1.797    0.085 .  
## phhown         0.735616    0.033766   10.450514  21.786 4.76e-10 ***
## phmov          0.034887    0.032953 7789.129217   1.059    0.290    
## phisp         -0.340442    0.019585 6914.823948 -17.383  < 2e-16 ***
## lhincm         0.521926    0.063343   85.086665   8.240 1.85e-12 ***
## punempm        0.394900    0.548156   80.339043   0.720    0.473    
## hwdis          0.357260    0.068071   75.061420   5.248 1.38e-06 ***
## fbdis         -0.137827    0.288961   82.215359  -0.477    0.635    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## singular fit

And test whether variation in the slopes and intercept is significant from 0.

ranova(re.int.slope.model)

## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## lhinc ~ concd + pfb + phhown + phmov + phisp + lhincm + punempm + 
##     hwdis + fbdis + (1 + phhown + pfb | GEOIDm)
##                                       npar  logLik    AIC    LRT Df
## <none>                                  17 -1843.7 3721.5          
## phhown in (1 + phhown + pfb | GEOIDm)   14 -1853.1 3734.2 18.685  3
## pfb in (1 + phhown + pfb | GEOIDm)      14 -1887.8 3803.6 88.161  3
##                                       Pr(>Chisq)    
## <none>                                              
## phhown in (1 + phhown + pfb | GEOIDm)  0.0003176 ***
## pfb in (1 + phhown + pfb | GEOIDm)     < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We can add random slopes to other variables. For example, we found that percent foreign born is positively associated with log Hispanic median household income. Does that vary across metro areas?

re.int.slope.model2 <- lmer(lhinc~ concd  +pfb + phhown   +phmov  + phisp 
              +lhincm + punempm + hwdis + fbdis + (1+phhown+pfb|GEOIDm), 
              data=ca.metro.tracts)
ranova(re.int.slope.model2)

## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## lhinc ~ concd + pfb + phhown + phmov + phisp + lhincm + punempm + 
##     hwdis + fbdis + (1 + phhown + pfb | GEOIDm)
##                                       npar  logLik    AIC    LRT Df
## <none>                                  17 -1843.7 3721.5          
## phhown in (1 + phhown + pfb | GEOIDm)   14 -1853.1 3734.2 18.685  3
## pfb in (1 + phhown + pfb | GEOIDm)      14 -1887.8 3803.6 88.161  3
##                                       Pr(>Chisq)    
## <none>                                              
## phhown in (1 + phhown + pfb | GEOIDm)  0.0003176 ***
## pfb in (1 + phhown + pfb | GEOIDm)     < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Visualizing differences

We can visualize the differences between the single level, random intercept, random coefficient, and random intercept and coefficient models by extracting the predicted values of log Hispanic median household income from each model. This will help fortify how exactly these different approaches are modelling the outcome. For simplicity’s sake, we’ll fit models only examining the relationship between percent Hispanic home ownership and log Hispanic median income. First, fit the models

ols <- lm(lhinc ~ phhown, data = ca.metro.tracts)
reint <- lmer(lhinc~ phhown  + (1|GEOIDm), data=ca.metro.tracts)
reslope <- lmer(lhinc~ (0+phhown|GEOIDm), data=ca.metro.tracts)
reintslope <- lmer(lhinc~  (1+phhown|GEOIDm), data=ca.metro.tracts)

Second, we extract the predicted or fitted log Hispanic median household income using the fitted() function. Let’s save these predictions into our data set ca.metro.tracts

ca.metro.tracts$PooledPredictions <- fitted(ols)
ca.metro.tracts$VaryingInterceptPredictions <- fitted(reint)
ca.metro.tracts$VaryingSlopePredictions <- fitted(reslope)
ca.metro.tracts$InteractionPredictions <- fitted(reintslope)

Next, we plot the random intercept predictions compared to the single level model predictions

ggplot(ca.metro.tracts) +
    geom_line(aes(x = phhown, y = PooledPredictions), color = "darkgrey") +
    geom_line(aes(x = phhown, y = VaryingInterceptPredictions), color = "blue") +
    #geom_line(aes(x = phhown, y = VaryingSlopePredictions), color = "red") +
    #geom_line(aes(x = phhown, y = InteractionPredictions), color = "black") +
    facet_wrap(~GEOIDm) +
    theme_bw()

Do the plots jibe with what you thought these models are doing?

Next, varying slope for phhown vs. single level

ggplot(ca.metro.tracts) +
    geom_line(aes(x = phhown, y = PooledPredictions), color = "darkgrey") +
    #geom_line(aes(x = phhown, y = VaryingInterceptPredictions), color = "blue") +
    geom_line(aes(x = phhown, y = VaryingSlopePredictions), color = "red") +
    #geom_line(aes(x = phhown, y = InteractionPredictions), color = "black") +
    facet_wrap(~GEOIDm) +
    theme_bw()

Finally, the varying slope and intercept vs single level

ggplot(ca.metro.tracts) +
    geom_line(aes(x = phhown, y = PooledPredictions), color = "darkgrey") +
    #geom_line(aes(x = phhown, y = VaryingInterceptPredictions), color = "blue") +
    #geom_line(aes(x = phhown, y = VaryingSlopePredictions), color = "red") +
    geom_line(aes(x = phhown, y = InteractionPredictions), color = "black") +
    facet_wrap(~GEOIDm) +
    theme_bw()

Based on the random slope + intercept model, which metro areas exhibit a phhown slope that considerably deviates from the mean? Which metro areas do not?

Single level vs. Multilevel

The ranova() function allows us to compare multilevel models with single level regressions. We already used ranova() above to test whether the variance of the random slope is statistically different from a null of 0. This test is the same thing as comparing model fit for a random slope model compared to a regular single level regression. We test the random intercept model against the single level as follows.

ranova(re.int.model3)

## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## lhinc ~ concd + pfb + phhown + phmov + phisp + lhincm + punempm + 
##     hwdis + fbdis + (1 | GEOIDm)
##              npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>         12 -1896.5 3817.1                         
## (1 | GEOIDm)   11 -1922.2 3866.5 51.377  1  7.623e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

What about the random slope + intercept model?

ranova(re.int.slope.model)

## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## lhinc ~ concd + pfb + phhown + phmov + phisp + lhincm + punempm + 
##     hwdis + fbdis + (1 + phhown + pfb | GEOIDm)
##                                       npar  logLik    AIC    LRT Df
## <none>                                  17 -1843.7 3721.5          
## phhown in (1 + phhown + pfb | GEOIDm)   14 -1853.1 3734.2 18.685  3
## pfb in (1 + phhown + pfb | GEOIDm)      14 -1887.8 3803.6 88.161  3
##                                       Pr(>Chisq)    
## <none>                                              
## phhown in (1 + phhown + pfb | GEOIDm)  0.0003176 ***
## pfb in (1 + phhown + pfb | GEOIDm)     < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The results actually compare models with random slopes for phhown and pfb separately. To test the combined model against a single level model, use the reduce.terms = FALSE argument

ranova(re.int.slope.model, reduce.terms = FALSE)

## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## lhinc ~ concd + pfb + phhown + phmov + phisp + lhincm + punempm + 
##     hwdis + fbdis + (1 + phhown + pfb | GEOIDm)
##                             npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                        17 -1843.7 3721.5                         
## (1 + phhown + pfb | GEOIDm)   11 -1922.2 3866.5 157.01  6  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Random intercept vs. Random slope

Random intercept and slope models are nested, which means we can use a basic ANOVA to make comparisons between the two. We do this by using the anova() function, which we also used when we made model comparisons between spatial regression models in Week 7 Lab.

The null will always be the simpler model, with a random intercept only model being the simplest among the multilevel flavors. For example, testing model fit between a random intercept only model and a random slope on phhown yields :

anova(re.int.model3, re.slope.model1)

## Data: ca.metro.tracts
## Models:
## re.int.model3: lhinc ~ concd + pfb + phhown + phmov + phisp + lhincm + punempm + 
## re.int.model3:     hwdis + fbdis + (1 | GEOIDm)
## re.slope.model1: lhinc ~ concd + pfb + phhown + phmov + phisp + lhincm + punempm + 
## re.slope.model1:     hwdis + fbdis + (0 + phhown | GEOIDm)
##                 Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## re.int.model3   12 3772.3 3855.9 -1874.2   3748.3                        
## re.slope.model1 12 3788.1 3871.7 -1882.1   3764.1     0      0          1

Random intercept only vs Random intercept and random slopes for phhown and pfb?

anova(re.int.model3, re.int.slope.model)

## Data: ca.metro.tracts
## Models:
## re.int.model3: lhinc ~ concd + pfb + phhown + phmov + phisp + lhincm + punempm + 
## re.int.model3:     hwdis + fbdis + (1 | GEOIDm)
## re.int.slope.model: lhinc ~ concd + pfb + phhown + phmov + phisp + lhincm + punempm + 
## re.int.slope.model:     hwdis + fbdis + (1 + phhown + pfb | GEOIDm)
##                    Df    AIC    BIC  logLik deviance  Chisq Chi Df
## re.int.model3      12 3772.3 3855.9 -1874.2   3748.3              
## re.int.slope.model 17 3675.0 3793.4 -1820.5   3641.0 107.34      5
##                    Pr(>Chisq)    
## re.int.model3                    
## re.int.slope.model  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The difference between ranova() and anova() is that the latter generalizes what the former is doing. ranova() tests against the null of a single level model. anova() allows you test across all nested models.

We can also rely on our new best friend the Akaike Information Criterion (AIC). Remember, AIC is related to deviance. Deviance is bad. So, the lower the deviance, the better, which means the lower the AIC, the better the fit. We can use the AIC to even compare across fixed and random effects models.

#Save AIC values
AICs<-c(AIC(single.model2),AIC(fe.int.model1), AIC(re.int.model3), AIC(re.slope.model1), AIC(re.int.slope.model))

#plot the AICs
plot(AICs, type="l", lwd=1.5, xaxt="n", xlab="")
axis(1, at=1:5,labels=F) #6= number of models
labels<-c("OLS", "FE","RE-Int", "RE-Slope","RE-Int/Slope" )
text(1:5, par("usr")[3]-.25, srt=45, adj=1, labels=labels, xpd=T)
mtext(side=1, text="Model Specification", line=3)
#circle the model with the lowest AIC
symbols(x= which.min(AICs), y=AICs[which.min(AICs)], circles=1, fg=2,lwd=2,add=T)

The random intercept and slope model seems to be the best fitting. But, once again, what were trying to do with an FE model is very different with what we are trying to accomplish with a multilevel model. What is your hypothesis, what is your theoretical framework, what pathways are you interested in…